Kim ANTUNEZ - Julien PRAMIL
octobre 2022
Plan :
btb (PSAR analyse urbaine, Arlindo Dos Santo et François Sémécurbe).kernelSmoothing, permet de réaliser très facilement un carroyage et un lissage sur des données géolocalisées avec R.À partir de données ponctuelles, nous allons apprendre en utilisant le langage R :
Liens utiles
Avant toute diffusion auprès des partenaires, il faut bien veiller à respecter :
Pour simplifier : on prend des libertés avec la sémiologie cartographique
Auteur : Timothée Giraud, auteur de la librairie mapsf
| Nom | Description | Code EPSG |
|---|---|---|
| Lambert93 | Système de projection officiel pour la métropole | 2154 |
| LAEA | Système de projection européen | 3035 |
| WGS84 | GPS (utile pour utiliser Leaflet) | 4326 |
Cinq librairies sont nécessaires pour ce TP.
sf pour manipuler des fichiers spatiaux (importer des .shp, transformer des projections, et réaliser des géotraitements) ;
mapsf pour réaliser des cartes dans RStudio ;
mapview (reposant sur leaflet) pour réaliser des cartes interactives (fond de carte OpenStreetMap);
btb pour le carroyage et lissage ;
dplyr pour le traitement des données, en particulier l’agrégation géographique.
Charger les librairies nécessaires
## Liste des librairies utilisées
packages <- c("dplyr","sf","btb","mapsf","leaflet","mapview")
## Vérifier si la librairie est installée, si non l'installer, puis la charger
package.check <- lapply(
packages,
FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE, quiet = TRUE)
library(x, character.only = TRUE)
}
}
)Sur le SSPCloud :
aws.s3 pour charger les données stockées dans Minio ;Elle recense :
À partir de cette source, nous avons constitué une base :
Avec les 8 variables suivantes :
id_mutation : identifiant unique de la ventedate_mutation : date de la ventetype_local : appartement ou maisonnombre_pieces_principales : nombre de pièces dans le logementvaleur_fonciere : prix de ventesurface_reelle_bati : surface du logementx : longitude (en projection Lambert 93)y : latitude (en projection Lambert 93)Remarques :
Base non représentative de la réalité : usage purement pédagogique.
Dans le SSPcloud :
Importation de la base ventesImmo_couronneParis.RDS, stockée sous Minio
En dehors du SSPCloud : téléchargement via URL.
On peut ensuite manipuler notre base de données chargée en mémoire.
id_mutation date_mutation type_local nombre_pieces_principales
1 2021-447023 2021-01-08 Appartement 3
2 2021-447024 2021-01-05 Appartement 2
3 2021-447025 2021-01-08 Appartement 3
4 2021-447027 2021-01-08 Appartement 3
5 2021-447029 2021-01-05 Appartement 2
6 2021-447030 2021-01-15 Appartement 3
valeur_fonciere surface_reelle_bati x y
1 480000 64 647357.3 6868635
2 345000 43 644483.5 6867695
3 384000 41 648001.8 6866153
4 261900 44 647414.8 6868937
5 407200 24 646929.9 6864730
6 1040000 90 645132.7 6864646
[1] 34489 8
.gpkg ou .shp
.gpkg (voir détails ici)Récupérons donc le contour du territoire étudié avec la fonction st_read du package sf.
chemin_file <- paste0(url_bucket,bucket,"/r-lissage-spatial/depCouronne.gpkg")
depCouronne_sf <- sf::st_read(chemin_file)Reading layer `depCouronne' from data source
`https://minio.lab.sspcloud.fr/projet-formation/r-lissage-spatial/depCouronne.gpkg'
using driver `GPKG'
Simple feature collection with 4 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
Projected CRS: RGF93 / Lambert-93
Visualisons cette couche vectorielle
Simple feature collection with 4 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
Projected CRS: RGF93 / Lambert-93
code libelle reg surf geom
1 75 Paris 11 105 MULTIPOLYGON (((660897 6860...
2 92 Hauts-de-Seine 11 176 MULTIPOLYGON (((648796 6847...
3 93 Seine-Saint-Denis 11 237 MULTIPOLYGON (((659428 6861...
4 94 Val-de-Marne 11 245 MULTIPOLYGON (((656908 6846...
On renomme la variable geom en geometry
Une fois les départements chargés, nous allons sélectionner le contour de la commune de Paris.
Visualisons cette couche vectorielle
Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 643076 ymin: 6857499 xmax: 660897 ymax: 6867034
Projected CRS: RGF93 / Lambert-93
code libelle reg surf geometry
1 75 Paris 11 105 MULTIPOLYGON (((660897 6860...
On reprojète dans le même système de projection (si besoin). Ici, en Lambert 93 (epsg 2154)
Pour éviter les effets de bord ==> sélectionner des données au-delà de notre zone d’intérêt (ici Paris intramuros)
Création d’un buffer de la bbox, avec une marge de 2000 mètres
Schéma explicatif des zones sélectionnées :
# Petit rectangle vectoriel
bbox_sf = st_sf(nom="bbox", geometry = st_as_sfc(bbox), crs = 2154)
# Grand rectangle vectoriel
bufferBbox_sf = st_sf(nom="buffer_bbox", geometry = st_as_sfc(bufferBbox), crs = 2154)
# Options globales pour les cartes avec mapview
mapviewOptions(
basemaps = c(
"OpenStreetMap",
"CartoDB.Positron",
"CartoDB.DarkMatter",
"Esri.WorldImagery"
)
)
# Cartographie pédagogique avec mapview
m <- mapview(paris_sf ,col.regions= "#26cce7")+
mapview(bufferBbox_sf %>% st_cast("MULTILINESTRING"),color="#FFC300",lwd=6)+
mapview(bbox_sf %>% st_cast("MULTILINESTRING"),color="#229954",lwd=6)
mFiltrer (numériquement) les logements dans le triangle jaune
xmin ymin xmax ymax
641076 6855499 662897 6869034
xMin <- bufferBbox["xmin"]
xMax <- bufferBbox["xmax"]
yMin <- bufferBbox["ymin"]
yMax <- bufferBbox["ymax"]
# Ne garder que les données dans le rectangle englobant, sans traitement vectoriel !
dfBase_filtre <- dfBase[dfBase$x >= xMin & dfBase$x <= xMax & dfBase$y >= yMin & dfBase$y <= yMax, ]
dim(dfBase_filtre)[1] 24419 8
[1] 34489 8
sf vectoriel ;Remarque: Pour la zone tampon, prendre une marge légèrement plus grande que le rayon de lissage envisagé.
Simple feature collection with 22221 features and 7 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 641099.2 ymin: 6855538 xmax: 662814.2 ymax: 6868990
Projected CRS: RGF93 / Lambert-93
First 10 features:
id_mutation date_mutation type_local nombre_pieces_principales
3 2021-447025 2021-01-08 Appartement 3
5 2021-447029 2021-01-05 Appartement 2
6 2021-447030 2021-01-15 Appartement 3
9 2021-447034 2021-01-18 Appartement 4
10 2021-447035 2021-01-18 Appartement 4
11 2021-447037 2021-01-07 Appartement 4
14 2021-447040 2021-01-07 Appartement 4
15 2021-447041 2021-01-04 Appartement 4
16 2021-447042 2021-01-15 Appartement 3
19 2021-447046 2021-01-14 Appartement 1
valeur_fonciere surface_reelle_bati nom geometry
3 384000 41 territoire POINT (648001.8 6866153)
5 407200 24 territoire POINT (646929.9 6864730)
6 1040000 90 territoire POINT (645132.7 6864646)
9 625000 68 territoire POINT (647494.6 6866019)
10 721250 94 territoire POINT (643891.9 6864587)
11 1035000 94 territoire POINT (648239 6866235)
14 652000 62 territoire POINT (647837.9 6865935)
15 830000 94 territoire POINT (643093.5 6864422)
16 407070 55 territoire POINT (647482.2 6867646)
19 235000 22 territoire POINT (647066.8 6866298)
Ci-dessous le code permettant d’avoir un aperçu de la zone, du buffer et de 2000 points tirés aléatoirement dedans.
# Mise en forme de la couche buffer
buffer_sf$nom <- "buffer"
# Échantillon de 2000 observations dans le buffer
sfBase_sample <- sfBase_filtre[sample(1:nrow(sfBase_filtre),2000) ,]
# Cartographie pédagogique avec mapview
mapview(paris_sf ,col.regions= "#26cce7")+
mapview(buffer_sf %>% st_cast("MULTILINESTRING"),color="#FFC300",lwd=6)+
mapview(sfBase_sample,#col.regions = "black",alpha.regions=0.5,
alpha=0,cex=2)Avant de lisser les données ponctuelles, on peut souhaiter représenter ces données sous forme carroyée afin de se les approprier. Il faut pour cela :
iCellSize = 200 # Carreaux de 200m
points_carroyage <- dfBase_filtre # On repart de la base filtrée selon la première méthode
points_carroyage$x_centroide = points_carroyage$x -
(points_carroyage$x %% iCellSize) + (iCellSize / 2)
points_carroyage$y_centroide = points_carroyage$y -
(points_carroyage$y %% iCellSize) + (iCellSize / 2)
head(points_carroyage) id_mutation date_mutation type_local nombre_pieces_principales
1 2021-447023 2021-01-08 Appartement 3
2 2021-447024 2021-01-05 Appartement 2
3 2021-447025 2021-01-08 Appartement 3
4 2021-447027 2021-01-08 Appartement 3
5 2021-447029 2021-01-05 Appartement 2
6 2021-447030 2021-01-15 Appartement 3
valeur_fonciere surface_reelle_bati x y x_centroide y_centroide
1 480000 64 647357.3 6868635 647300 6868700
2 345000 43 644483.5 6867695 644500 6867700
3 384000 41 648001.8 6866153 648100 6866100
4 261900 44 647414.8 6868937 647500 6868900
5 407200 24 646929.9 6864730 646900 6864700
6 1040000 90 645132.7 6864646 645100 6864700
Remarque : Ces deux premières étapes ne sont pas intégrées dans des fonctions du package btb. Il faut donc s’inspirer du code ci-dessus pour reproduire le même type de carroyage. Il est envisagé qu’une future version du package btb simplifie le code permettant de réaliser un carroyage.
Utilisation de dfToGrid du package btb, avec comme paramètres obligatoires.
df : un tableau avec les colonnes x et y représentant les coordonnées des centroïdes de la grille ;sEPSG : une chaîne de caractères indiquant le code epsg du système de projection utilisé ;iCellSize : la taille des carreaux (longueur du côté, en mètres).On obtient le carroyage des ventes dans Paris intramuros
contourParis <- st_cast(paris_sf[,c("geometry")],"MULTILINESTRING")
mf_init(x=carreaux,theme = "agolalight")
mf_map(x = carreaux,
type = "choro",
var="nbVentes",
breaks = "quantile",
nbreaks = 5,
lwd=1,
leg_val_rnd = 1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Carroyage du nombre de ventes",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Remarque :
Le carroyage pourrait aussi être utilisé pour simplifier les données avant lissage si le calcul de lissage est trop long en raison d’un très grand nombre d’observations.
Toujours pour éviter les effets de bord :
dfBase_filtre (construite précédemment).nbObsLisse qui vaudra 1 pour chaque observationnbObsLisse avec la fonction kernelSmoothing du package btb. id_mutation date_mutation type_local nombre_pieces_principales
1 2021-447023 2021-01-08 Appartement 3
2 2021-447024 2021-01-05 Appartement 2
3 2021-447025 2021-01-08 Appartement 3
4 2021-447027 2021-01-08 Appartement 3
5 2021-447029 2021-01-05 Appartement 2
6 2021-447030 2021-01-15 Appartement 3
valeur_fonciere surface_reelle_bati x y nbObsLisse
1 480000 64 647357.3 6868635 1
2 345000 43 644483.5 6867695 1
3 384000 41 648001.8 6866153 1
4 261900 44 647414.8 6868937 1
5 407200 24 646929.9 6864730 1
6 1040000 90 645132.7 6864646 1
Paramètres de kernelSmoothing :
dfObservations : le tableau des données à lisser. Il doit nécessairement contenir une colonne x, une colonne y, et 1 à n colonnes numériques (variables à lisser) ;sEPSG : chaine de caractères indiquant le code epsg du système de projection utilisé ;iCellSize : un entier indiquant la taille des carreaux ;iBandwidth : un entier indiquant le rayon de lissage.Lissage et enregistrement dans sfCarrLiss
Visualisation du résultat :
Simple feature collection with 6 features and 3 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 646000 ymin: 6855400 xmax: 647800 ymax: 6855600
Projected CRS: RGF93 / Lambert-93
x y nbObsLisse geometry
1 646100 6855500 1.712851 POLYGON ((646000 6855600, 6...
2 646500 6855500 1.261316 POLYGON ((646400 6855600, 6...
3 647100 6855500 2.301557 POLYGON ((647000 6855600, 6...
4 647300 6855500 2.808277 POLYGON ((647200 6855600, 6...
5 647500 6855500 2.395913 POLYGON ((647400 6855600, 6...
6 647700 6855500 1.949819 POLYGON ((647600 6855600, 6...
Nombre de lignes lissées
Remarque :
Attention, la fonction retourne une erreur en cas de :
x ou y.La fonction de lissage classique du package btb est conservative.
Remarque sur le secret :
On cartographie les carreaux :
choroplèthe ;# Filtrage des carreaux lissés intersectant la ville de Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf[,"geometry"],left=F)
# Carte lissée
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 400m",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")On souhaite obtenir le nombre de ventes lissé, en 2021 dans Paris intramuros :
kernelSmoothing
nbObsLissee)Taille optimale du rayon de lissage ?
Inconnu a priori…
Compromis entre précision et généralisation
On teste avec :
# Carte lissée
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 600m",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")# Carte lissée
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Le choix revient donc au statisticien-géographe, en fonction de sa connaissance des données et des objectifs recherchés.
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
border = NA, # C'est ici que ça se passe
# lwd=1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m, sans grille",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Arètes des carreaux visibles
Idée : réduite la taille des carreaux (=> 50m)
Attention au temps de calcul !
(et on en profite aussi pour faire varier le rayon de lissage => 1km)
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
border = NA,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec pas de 50m",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Que se passe-t-il si on lisse uniquement les ventes réalisées dans Paris intramuros ?
Comparons avec les mêmes paramètres :
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
type = "choro",
var="nbObsLisse",
breaks = c(0,3,5,7,9,17),
# nbreaks = 5,
border = NA,
leg_val_rnd = 1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Sans effets de bord, avec R=2000m",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")sfBase_intramuros <- sfBase %>% st_join(paris_sf,left=F)
dfBase_intramuros <- sfBase_intramuros %>%
mutate(x=st_coordinates(geometry)[,1],
y=st_coordinates(geometry)[,2],
nbObsLisse=1) %>%
st_drop_geometry() %>%
select(nbObsLisse,x,y)
sfCarLiss_intramuros <- btb::kernelSmoothing(dfObservations = dfBase_intramuros,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 2000)
sfCarLiss_intramuros_paris <- sfCarLiss_intramuros %>%
st_join(paris_sf,left = F)mf_init(x=sfCarLiss_intramuros_paris,theme = "agolalight")
mf_map(x = sfCarLiss_intramuros_paris,
type = "choro",
var="nbObsLisse",
breaks = c(0,3,5,7,9,17),
# nbreaks = 5,
border = NA,
leg_val_rnd = 1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Avec effets de bord, avec R=2000m",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Les effets de bord se matérialisent à la périphérie de Paris, où les densités lissées sont artificiellement plus faibles que quand on lisse en prenant une marge.
btb
Par défaut : la grille produite automatiquement par btb::kernelSmoothing dépasse les limites de la zone d’étude choisie.
Prenons l’exemple précédent du lissage sur les ventes à l’intérieur de Paris intramuros.
btb
La grille obtenue dépasse Paris intramuros et contient :
btb
Modification de cette règle avec le paramètre iNeighbor.
Par exemple :
iNeighbor=2 => 2 carreaux limitrophes aux carreaux contenant au moins une vente.iNeighbor=0 => grille uniquement constituée des carreaux contenant au moins une vente
btb
btb
On peut souhaiter utiliser une grille de lissage particulière et adaptée à un territoire donné.
Par exemple, sur une ville cotière, on ne veut :
iNeighbor=0.Dans ce cas :
dfCentroids de la fonction kernelSmoothing.Dans cette section, nous allons nous intéresser au prix moyen des logements vendus à Paris en 2021.
Remarque :
Une moyenne n’est pas le lissage du rapport liss(variable / nbObsLisse) mais le rapport des lissages liss(variable) / liss(nbObs).
Il faut :
# Calculer le taux
sfCarLiss$prixMoyen <- sfCarLiss$valeur_fonciere / sfCarLiss$nbObsLisse
# Cartographie
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf,left=F)
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
type = "choro",
var="prixMoyen",
breaks = "quantile",
nbreaks = 5,
border=NA,
leg_val_rnd = 1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage des prix de vente avec un rayon de 800m",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Les prix moyens sont particulièrement élevés dans le centre et l’Ouest de la capitale
Remarque :
le lissage en moyenne est potentiellement sensible aux valeurs atypiques (contrairement au lissage quantile)
Nous nous intéressons ici à la proportion de ventes portant sur des grands logements (comportant plus de 4 pièces principales).
quatrePieces et nbObsLisse
# Lissage
dataLissage <- dfBase_filtre[,c("nbObsLisse","quatrePieces","x","y")]
sfCarLiss <- btb::kernelSmoothing(dfObservations = dataLissage,
sEPSG = "2154",
iCellSize = 100,
iBandwidth = 800)
# Calculer le taux
sfCarLiss$txQuatrePieces <- sfCarLiss$quatrePieces / sfCarLiss$nbObsLisse
# Intersection géographique avec Paris
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf,left=F)mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
type = "choro",
var="txQuatrePieces",
breaks = "quantile",
nbreaks = 5,
border=NA,
leg_val_rnd = 1,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Logements avec 4 pièces ou plus (rayon de 800m)",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Fort logiquement, la carte des prix moyens et la carte des grands logements ont des similitudes importantes.
On peut alors avoir envie de lisser le prix au mètre-carré.
Il convient de lisser séparément les prix de vente (numérateur) et le nombre de mètres-carrés (dénominateur) des logements vendus.
Autrement, si on lisse le prix au m² de chaque logement vendu, on sur-pondère artificiellement les prix au m² des petits logements (car l’unité statistique devient alors le logement, et non le m² vendu).
sfCarLiss_paris <- sfCarLiss[unlist(st_intersects(paris_sf,sfCarLiss)),]
mf_init(x = sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
type = "choro",
var="prixM2",
breaks = "quantile",
nbreaks = 5,
border=NA,
leg_val_rnd = 1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Prix au m² (rayon de 800m)",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")kernelSmoothing
vQuantiles, qui contient la liste des quantiles souhaitésContrairement au lissage classique :
nbObs indiquant le nombre réel (non lissé) d’observations ayant contribué au calcul du carreau.Nous allons calculer le 1er décile, la médiane et le 9ème décile du prix de vente des logements à Paris.
Visualisons les première lignes de la grille lissée
Simple feature collection with 6 features and 6 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 651500 ymin: 6854700 xmax: 652100 ymax: 6854800
Projected CRS: RGF93 / Lambert-93
nbObs valeur_fonciere_01 valeur_fonciere_05 valeur_fonciere_09 x y
1 57 182000 310000 610000 651550 6854750
2 53 182000 310000 620000 651650 6854750
3 52 182000 310000 620000 651750 6854750
4 55 182000 321000 620000 651850 6854750
5 48 182000 324000 620000 651950 6854750
6 55 182000 324000 620000 652050 6854750
geometry
1 POLYGON ((651500 6854800, 6...
2 POLYGON ((651600 6854800, 6...
3 POLYGON ((651700 6854800, 6...
4 POLYGON ((651800 6854800, 6...
5 POLYGON ((651900 6854800, 6...
6 POLYGON ((652000 6854800, 6...
Cartographie du 1er décile de prix des ventes
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf,left = F)
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
type = "choro",
var="valeur_fonciere_01",
breaks = "quantile",
nbreaks = 5,
border=NA,
leg_val_rnd = 0,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Premier décile des prix de vente (rayon de 1500m)",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Cartographie de la médiane des prix des ventes
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
type = "choro",
var="valeur_fonciere_05",
breaks = "quantile",
nbreaks = 5,
border=NA,
leg_val_rnd = 0,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Médiane des prix de vente (rayon de 1500m)",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Cartographie du rapport interdécile des prix des ventes
sfCarLiss_paris$rapp_interdecile <- sfCarLiss_paris$valeur_fonciere_09 / sfCarLiss_paris$valeur_fonciere_01
# Cartographie
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
type = "choro",
var="rapp_interdecile",
breaks = "quantile",
nbreaks = 5,
border=NA,
leg_val_rnd = 0,
add = TRUE)
mf_map(x = contourParis,
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Rapport interdéciles des prix de vente (rayon de 1500m)",
credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")Enfin, pour que le quantile lissé ait du sens, il faut qu’un nombre conséquent d’observations ait participé à sa construction. Ainsi, il est conseillé d’utiliser un rayon de lissage suffisamment élevé.
Objectif : agréger des données sur un zonage particulier
2 cas :
Calculons le nombre et le prix moyen des transactions immobilières en 2021 dans le triangle d’or.
Pour information, ce triangle a été créé manuellement sur le site du Géoportail.
.kml et se lit sans problème avec la fonction sf::st_read.On visualise le triangle :
[1] "WGS 84"
triangleOr <- st_transform(triangleOr,crs = 2154) # Transformation en Lambert93
# Intersection géographique avec les transactions
transac_triangle <- sfBase %>% st_join(triangleOr[,"geometry"],left=F)
head(transac_triangle,4) # Visualisation de points retenusSimple feature collection with 4 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 648814.5 ymin: 6863055 xmax: 649194 ymax: 6863601
Projected CRS: RGF93 / Lambert-93
id_mutation date_mutation type_local nombre_pieces_principales
18629 2021-489256 2021-01-15 Appartement 5
18676 2021-489336 2021-01-25 Appartement 3
18677 2021-489337 2021-01-28 Appartement 4
18698 2021-489369 2021-02-05 Appartement 6
valeur_fonciere surface_reelle_bati geometry
18629 3300000 180 POINT (648918.9 6863145)
18676 590000 53 POINT (648903.5 6863601)
18677 1511345 65 POINT (649194 6863377)
18698 2450000 133 POINT (648814.5 6863055)
Cartographie des points dans le triangle
Objectif : obtenir la proportion de logements sociaux dans un rayon de 1km autour de la Direction générale de l’Insee.
On utilise les données carroyées à 200 mètres disponibles sur insee.fr.
Remarque : un programme généralisant cet exercice est disponible sur insee.fr, en accompagnement des données carroyées.
# Fonction permettant de faire le chargement souhaité
st_read_maison <- function(chemin_tab){
requete <- "SELECT IdINSPIRE,Depcom,I_est_cr,Men, Log_soc, geom
FROM Filosofi2015_carreaux_200m_metropole
WHERE SUBSTR(Depcom, 1, 2) IN ('75','92','93','94') "
sf::st_read(chemin_tab, query = requete)
}
# Chargement des données
chemin_file <- paste0(url_bucket, bucket,
"/r-lissage-spatial/Filosofi2015_carreaux_200m_metropole.gpkg")
carreaux <- st_read_maison(chemin_file)
head(carreaux)Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 652623.2 ymin: 6858353 xmax: 655268.9 ymax: 6866380
Projected CRS: RGF93 / Lambert-93
IdINSPIRE Depcom I_est_cr Men Log_soc
1 CRS3035RES200mN2893400E3763200 75119 0 990 937
2 CRS3035RES200mN2890000E3762200 75111 0 926 80
3 CRS3035RES200mN2893400E3762400 75119 0 508 323
4 CRS3035RES200mN2891800E3763600 75119 0 633 0
5 CRS3035RES200mN2890800E3763200 75120 0 349 188
6 CRS3035RES200mN2885800E3760600 75113 0 751 344
geom
1 MULTIPOLYGON (((654521.9 68...
2 MULTIPOLYGON (((653843.4 68...
3 MULTIPOLYGON (((653725.1 68...
4 MULTIPOLYGON (((655069.7 68...
5 MULTIPOLYGON (((654764.7 68...
6 MULTIPOLYGON (((652641.8 68...
Vous êtes désormais capables de carroyer et lisser toutes les données que vous aurez à votre disposition !
BTBpy
